
;; spawn,'find ./irdc_dist_model/. -name "irdc_morph_pvec_*.sav" | xargs -I {} mv -v {} ./irdc_dist_model/pvec_data/.'

;;TRIM_MORPH_PVEC

s = read_bgps_csv()
restore,'./local/bgps_velocity_struct.sav',/ver
restore,'./irdc_dist_model/bgps_rb3.sav',/ver

alpha = findgen(51)/20.+0.5
salpha = string(alpha,format="(F0.2)")
fn = './irdc_dist_model/pvec_data/irdc_morph_pvec_'+salpha+'.sav'
n = n_elements(fn)

;; Get galactic params
defsysv, '!MW', exists = exists
IF ~ exists THEN galactic_params 
d = dindgen(!MW.NBINS)*!MW.BINSIZE + !MW.BINSTART


mean_FWHM  = fltarr(n)
mean_ratio = fltarr(n)
mean_pchoo = fltarr(n)

ratios = 10.^(findgen(51)/25)
nrat = n_elements(ratios)

pchoo_arr = dblarr(n,nrat)
rhist_arr = dblarr(n,nrat)
;; mean_rgn   = fltarr(n)
;; mean_rgy   = fltarr(n)
;; mean_fgn   = fltarr(n)
;; mean_fgy   = fltarr(n)
;; frac_yes   = fltarr(n)

;; Build the KDIST DPDFs here, then just use...
message,'Building the KDIST DPDFs...',/inf
restore,fn[0]
ind  = where_array(pvec.cnum, s.cnum, nhas)
l    = s[ind].glon_peak
b    = s[ind].glat_peak
v    = velocity[ind].densegas
vstd = velocity[0].v_std
str  = s[ind]
dtan = !MW.R0 * cos(l*!dtor)/cos(b*!dtor)

kdist = dindgen(nhas,n_elements(d))
ph2   = dindgen(nhas,n_elements(d))
FOR k=0,nhas-1 DO BEGIN
   kdist[k,*] = kdist_spectrum(l[k],b[k],vstd,v[*,k])
   ph2[k,*]   = prob_lonlat(str[k])
ENDFOR



message,'Starting to loop over the alphas...',/inf
;; Loop over the alphas...
FOR i=0L,n-1 DO BEGIN
   IF FILE_TEST(fn[i]) THEN BEGIN
      
      undefine,pvec
      restore,fn[i]
      
      nhas = n_elements(pvec)
      
      FWHM     = fltarr(nhas)
      ratio    = fltarr(nhas)
      pchoo    = fltarr(nhas)
      ;; ngrs     = long(total(rb3.grs))
      ;; ratio_gy = fltarr(ngrs)
      ;; ratio_gn = fltarr(ngrs)
      ;; fwhm_gy  = fltarr(ngrs)
      ;; fwhm_gn  = fltarr(ngrs)
      
      gii = 0
      
      ;; Loop over the objects...
      FOR j=0L,nhas-1 DO BEGIN
         
         ;; Calculate FWHM
         yfit = mpfitpeak(d,pvec[j].invx,A,nterms=4)
         FWHM[j] = 2.355 * A[2] / 1.d3
         
         ;; Calculate Peak Probabilities on either side of d_tan
         totprob = pvec[j].invx * kdist[j,*] * ph2[j,*]
         totprob /= max(totprob,maxind)
         
         nind = WHERE(d LE dtan[j],nnear)
         find = WHERE(d GT dtan[j],nfar)
         ;;print,nnear,nfar,dtan[j]
         nmax = max(totprob[nind])
         fmax = max(totprob[find])
         
         totprob /= total(totprob)
         
         IF nmax GE fmax THEN BEGIN
            ratio[j] = nmax/fmax 
            pchoo[j] = total(totprob[nind])
         ENDIF ELSE BEGIN
            ratio[j] = fmax/nmax 
            pchoo[j] = total(totprob[find])
         ENDELSE
         
         IF 0 THEN BEGIN
            
            dmorph = d[maxind]
            
            ;; plot,d/1.d3,totprob,xtit='D [kpc]',ytit='Prob'
            ;; oplot,d/1.d3,kdist[j,*]/max(kdist[j,*]),color=cgColor('Red')
            ;; vline,dmorph/1.d3,thick=3
            ;; vline,rb3[j].grs_dist,thick=3,color=cgColor('Green')
            ;; vline,dtan[j]/1.d3,thick=3,color=cgColor('Blue')
            ;; wait,0.5
            IF abs(dmorph/1.d3 - rb3[j].grs_dist) LE 1.0 THEN BEGIN
               ratio_gy[gii] = ratio[j]
               fwhm_gy[gii]  = FWHM[j]
            ENDIF ELSE BEGIN
               ratio_gn[gii] = ratio[j]
               fwhm_gn[gii]  = FWHM[j]
            ENDELSE
            
            gii++
            
         ENDIF
      ENDFOR                       ;; End OBJECT loop
      
      plot,ratio,pchoo,psym=cgsymcat(16),symsize=0.5,/xlog,$
           xtit='Ratio',ytit='P!dcnoo!n'
      oplot,ratios,ratios*0.+0.5,psym=2,color=cgColor('Cyan')
      
      sorti = sort(ratio)
      ratio = ratio[sorti]
      pchoo = pchoo[sorti]
      
      pchoo = median(pchoo,21)
      pfunc = interpol(pchoo,ratio,ratios) < 1.0
      oplot,ratios,pfunc,color=cgColor('Deep Pink'),thick=2
      
      rcdf = total(histogram(alog10(ratio),min=0.,binsize=0.001,loc=rh_l),/cum)
      rcdf /= max(rcdf)
      rhist = interpol(rcdf,10.^(rh_l),ratios)
      
      oplot,ratios,rhist,psym=10
      pchoo_arr[i,*] = pfunc
      rhist_arr[i,*] = rhist
      
      mean_FWHM[i]  = median(FWHM)
      mean_ratio[i] = median(ratio)
      mean_pchoo[i] = median(pchoo)
      ;; mean_rgy[i]   = median(ratio_gy[where(ratio_gy NE 0, ny)])
      ;; mean_rgn[i]   = median(ratio_gn[where(ratio_gn NE 0, nn)])
      ;; mean_fgy[i]   = median(fwhm_gy[where(fwhm_gy NE 0)])
      ;; mean_fgn[i]   = median(fwhm_gn[where(fwhm_gn NE 0)])
      ;; frac_yes[i]   = float(ny) / float(ny + nn)
      
      print,'Alpha='+salpha[i]+'  Mean FWHM='+$
            string(mean_FWHM[i],format="(F0.1)")+' kpc  Mean Ratio='+$
            string(mean_ratio[i],format="(F0.1)")+'  Mean P_choo='+$
            string(mean_pchoo[i],format="(F0.3)")
      
   ENDIF
ENDFOR                          ;; End ALPHA loop

save,alpha,mean_FWHM,mean_ratio,mean_pchoo,pchoo_arr,ratios,rhist_arr,$
     filename='./irdc_dist_model/beta_cube.sav',/ver



;; myps,'./irdc_dist_model/analysis_plots/alpha_testing_plots.eps',$
;;      xsize=16,ysize=5,/cmyk
;; symalph = cgSymbol('beta')
;; multiplot,[3,1],xgap=0.04,/doxaxis,mtitle='DPDF '+cgSymbol_extra('propto')+' ('+$
;;           cgSymbol('chi')+'!u2!n)!u-'+symalph+'!n'

;; ;;========
;; ;; Panel 1
;; plot,alpha,mean_FWHM,xtit=symalph,ytit='Median FWHM [kpc]',thick=3,yr=[0,10]
;; oplot,alpha,mean_fgy,thick=3,color=cgColor('Green')
;; oplot,alpha,mean_fgn,thick=3,color=cgColor('Red')

;; al_legend,/top,/right,box=0,['Entire Sample','GRS YES match','GRS NO match'],$
;;           linestyle=0,thick=3,color=['Opposite','Green','Red'],charsize=0.8

;; multiplot,/doyaxis,/doxaxis

;; ;;========
;; ;; Panel 2
;; plot,alpha,mean_ratio,xtit=symalph,ytit='Median Ratio of Peak Probabilities',$
;;      thick=3,yr=[0,40.]

;; oplot,alpha,mean_rgy,thick=3,color=cgColor('Green')
;; oplot,alpha,mean_rgn,thick=3,color=cgColor('Red')

;; al_legend,/top,/left,box=0,['Entire Sample','GRS YES match','GRS NO match'],$
;;           linestyle=0,thick=3,color=['Opposite','Green','Red'],charsize=0.8

;; multiplot,/doyaxis,/doxaxis

;; ;;========
;; ;; Panel 3
;; plot,alpha,frac_yes,thick=3,yr=[0,1],xtit=symalph,$
;;      ytit='GRS YES Fraction'

;; myps,/done,/mp

END
